{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div style='background-image: url(\"../../../share/images/header.svg\") ; padding: 0px ; background-size: cover ; border-radius: 5px ; height: 250px'>\n",
    "    <div style=\"float: right ; margin: 50px ; padding: 20px ; background: rgba(255 , 255 , 255 , 0.7) ; width: 50% ; height: 150px\">\n",
    "        <div style=\"position: relative ; top: 50% ; transform: translatey(-50%)\">\n",
    "            <div style=\"font-size: xx-large ; font-weight: 900 ; color: rgba(0 , 0 , 0 , 0.8) ; line-height: 100%\">Computational Seismology</div>\n",
    "            <div style=\"font-size: large ; padding-top: 20px ; color: rgba(0 , 0 , 0 , 0.5)\"> Lamb's problem </div>\n",
    "        </div>\n",
    "    </div>\n",
    "</div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<p style=\"width:20%;float:right;padding-left:50px\">\n",
    "<img src=../../../share/images/book.jpg>\n",
    "<span style=\"font-size:smaller\">\n",
    "</span>\n",
    "</p>\n",
    "\n",
    "\n",
    "---\n",
    "\n",
    "This notebook is part of the supplementary material \n",
    "to [Computational Seismology: A Practical Introduction](https://global.oup.com/academic/product/computational-seismology-9780198717416?cc=de&lang=en&#), \n",
    "Oxford University Press, 2016.\n",
    "\n",
    "\n",
    "##### Authors:\n",
    "* David Vargas ([@dvargas](https://github.com/davofis))\n",
    "* Heiner Igel ([@heinerigel](https://github.com/heinerigel))\n",
    "\n",
    "\n",
    "(Based on the original code by Lane Johnson)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic Equations\n",
    "\n",
    "The fundamental analytical solution to the three-dimensional Lamb’s problem, the problem of determining the elastic disturbance resulting from a point force in a homogeneous half space, is implemented in this Ipython notebook. This solution provides fundamental information used as benchmark for comparison with entirely numerical solutions. A setup of the fundamental problem is illustrated below. The figure on the right hand side is published in [1] (Figure 1. System of coordinates)\n",
    "\n",
    "<p style=\"width:65%;float:right;padding-left:50px\">\n",
    "<img src=lambs_setup.png>\n",
    "<span style=\"font-size:smaller\">\n",
    "</span>\n",
    "</p>\n",
    "\n",
    "Simulations of 3D elastic wave propagation need to be validated by the use of analytical solutions. In order to evaluate how healthy a numerical solution is, one may recreate conditions for which analytical solutions exist with the aim of reproducing and compare the different results.\n",
    "\n",
    "We which to find the displacement wavefield $\\mathbf{u}(\\mathbf{x},t)$ at some distance $\\mathbf{x}$ from a seismic  source with $ \\mathbf{F} = f_1\\mathbf{\\hat{x}_1} + f_2\\mathbf{\\hat{x}_2} + f_3\\mathbf{\\hat{x}_3}$.\n",
    "\n",
    "For a uniform elastic material and a Cartesian co-ordinate system the equation for the conservation of linear momentum can be written\n",
    "\n",
    "\\begin{align*}\n",
    "\\rho(x) \\frac{\\partial^2}{\\partial t^2} \\mathbf{u(\\mathbf{x},t)} = (\\lambda + \\mu)\\nabla(\\nabla\\mathbf{u(\\mathbf{x},t)}) + \\mu\\nabla^2 \\mathbf{u(\\mathbf{x},t)} + \\mathbf{f(\\mathbf{x},t)}\n",
    "\\end{align*}\n",
    "\n",
    "We will consider the case where the source function is localized in both time and space\n",
    "\n",
    "\\begin{align*}\n",
    "\\mathbf{f(\\mathbf{x},t)} = (f_1\\mathbf{\\hat{x}_1} + f_2\\mathbf{\\hat{x}_2} + f_3\\mathbf{\\hat{x}_3})\\delta(x_1 - x^{'}_{1})\\delta(x_2 - x^{'}_{2})\\delta(x_3 - x^{'}_{3})\\delta(t - t^{'})\n",
    "\\end{align*}\n",
    "\n",
    "For such a source we will refer to the displacement solution as a Green’s function, and use the standard notation\n",
    "\n",
    "\\begin{align*}\n",
    "\\mathbf{u(\\mathbf{x},t)} = g_1(\\mathbf{x},t;\\mathbf{x^{'}},t^{'})\\mathbf{\\hat{x}_1} + g_2(\\mathbf{x},t;\\mathbf{x^{'}},t^{'})\\mathbf{\\hat{x}_2} + g_3(\\mathbf{x},t;\\mathbf{x^{'}},t^{'})\\mathbf{\\hat{x}_3}\n",
    "\\end{align*}\n",
    "\n",
    "The complete solution is found after applying the Laplace transform to the elastic wave equation, implementing the stress-free boundary condition, defining some transformations, and performing some algebraic manoeuvres. Then, the Green's function at the free surface is given:\n",
    "\n",
    "\\begin{align*}\n",
    "\\begin{split}\n",
    "\\mathbf{G}(x_1,x_2,0,t;0,0,x^{'}_{3},0) & =  \\dfrac{1}{\\pi^2\\mu r} \\dfrac{\\partial}{\\partial t}\\int_{0}^{((t/r)^2 - \\alpha^{-2})^{1/2}}\\mathbf{H}(t-r/\\alpha)\\mathbb{R}[\\eta_\\alpha\\sigma^{-1}((t/r)^2 - \\alpha^{-2} - p^2)^{-1/2}\\mathbf{M}(q,p,0,t,x^{'}_{3})\\mathbf{F}] dp \\\\\n",
    " & + \\dfrac{1}{\\pi^2\\mu r} \\dfrac{\\partial}{\\partial t}\\int_{0}^{p_2}\\mathbf{H}(t-t_2)\\mathbb{R}[\\eta_\\beta\\sigma^{-1}((t/r)^2 - \\beta^{-2} - p^2)^{-1/2}\\mathbf{N}(q,p,0,t,x^{'}_{3})\\mathbf{F}] dp\n",
    "\\end{split}\n",
    "\\end{align*}\n",
    "\n",
    "Details on the involved terms are found in the original paper [2]. The Green's $\\mathbf{G}$ function consist of three components of displacement evolving from the application of three components of force $\\mathbf{F}$. If we assume that each component of $\\mathbf{F}$ provokes three components of displacement, then $\\mathbf{G}$ is composed by nine independent components that correspond one to one to the matrices $\\mathbf{M}$ and $\\mathbf{N}$. Without losing generality it is shown that among them four are equal zero, and we end up only with five possible components. \n",
    "\n",
    "<p style=\"text-align: justify;\">\n",
    " [1] Eduardo Kausel - Lamb's problem at its simplest, 2012</p>\n",
    " \n",
    "<p style=\"text-align: justify;\">\n",
    " [2] Lane R. Johnson - Green’s Function for Lamb’s Problem, 1974</p>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Import all necessary libraries, this is a configuration step for the exercise.\n",
    "# Please run it before the simulation code!\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "from ricker import ricker\n",
    "\n",
    "# Show the plots in the Notebook.\n",
    "plt.switch_backend(\"nbagg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Compile the source code (needs gfortran!)\n",
    "!rm -rf lamb.exe output.txt\n",
    "!gfortran canhfs.for -fcheck=all -o lamb.exe"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calling the original FORTRAN code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Initialization of setup:\n",
    "# Figure 4 in Lane R. Johnson - Green’s Function for Lamb’s Problem, 1974\n",
    "# is reproduced when the following parameters are given\n",
    "# -----------------------------------------------------------------------------\n",
    "r   = 10.0    # km\n",
    "vp  = 8.0     # P-wave velocity km/s\n",
    "vs  = 4.62    # s-wave velocity km/s\n",
    "rho = 3.3     # Density kg/m^3 \n",
    "nt  = 512     # Number of time steps\n",
    "dt  = 0.01    # Time step s\n",
    "h   = 0.2     # Source position km (0.01 to reproduce Fig 2.16 of the book)\n",
    "ti  = 0.0     # Initial time s\n",
    "\n",
    "var = [vp, vs, rho, nt, dt, h, r, ti]\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# Execute fortran code\n",
    "# -----------------------------------------------------------------------------\n",
    "with open('input.txt', 'w') as f:\n",
    "    for i in var:\n",
    "        print(i, file=f, end='  ') # Write input for fortran code\n",
    "\n",
    "f.close()\n",
    "\n",
    "os.system(\"./lamb.exe\")   # Code execution\n",
    "\n",
    "# -----------------------------------------------------------------------------\n",
    "# Load the solution\n",
    "# -----------------------------------------------------------------------------\n",
    "G = np.genfromtxt('output.txt')\n",
    "\n",
    "u_rx = G[:,0]    # Radial displacement owing to horizontal load\n",
    "u_tx = G[:,1]    # Tangential displacement due to horizontal load\n",
    "u_zx = G[:,2]    # Vertical displacement owing to horizontal load\n",
    "\n",
    "u_rz = G[:,3]    # Radial displacement owing to a vertical load\n",
    "u_zz = G[:,4]    # Vertical displacement owing to vertical load\n",
    "\n",
    "t = np.linspace(dt, nt*dt, nt)    # Time axis\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualization of the Green's function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "code_folding": [
     0
    ],
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Plotting\n",
    "# -----------------------------------------------------------------------------\n",
    "seis = [u_rx, u_tx, u_zx, u_rz, u_zz]  # Collection of seismograms\n",
    "labels = ['$u_{rx}(t) [cm]$','$u_{tx}(t)[cm]$','$u_{zx}(t)[cm]$','$u_{rz}(t)[cm]$','$u_{zz}(t)[cm]$']\n",
    "cols = ['b','r','k','g','c']\n",
    "\n",
    "# Initialize animated plot\n",
    "fig = plt.figure(figsize=(12,8), dpi=80)\n",
    "\n",
    "fig.suptitle(\"Green's Function for Lamb's problem\", fontsize=16)\n",
    "\n",
    "plt.ion() # set interective mode\n",
    "plt.show()\n",
    "\n",
    "for i in range(5):              \n",
    "    st = seis[i]\n",
    "    ax = fig.add_subplot(2, 3, i+1)\n",
    "    ax.plot(t, st, lw = 1.5, color=cols[i])  \n",
    "    ax.set_xlabel('Time(s)')\n",
    "    ax.text(0.8*nt*dt, 0.8*max(st), labels[i], fontsize=16)\n",
    "    plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))\n",
    "    \n",
    "    ax.spines['left'].set_position('zero')\n",
    "    ax.spines['right'].set_color('none')\n",
    "    ax.spines['bottom'].set_position('zero')\n",
    "    ax.spines['top'].set_color('none')\n",
    "    ax.spines['left'].set_smart_bounds(True)\n",
    "    ax.spines['bottom'].set_smart_bounds(True)\n",
    "    ax.xaxis.set_ticks_position('bottom')\n",
    "    ax.yaxis.set_ticks_position('left')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Convolution\n",
    "\n",
    "Let $S(t)$ be a general source time function, then the displacent seismogram is given in terms of the Green's function $G$ via\n",
    "\n",
    "\\begin{equation}\n",
    "u(\\mathbf{x},t) = G(\\mathbf{x},t; \\mathbf{x}',t') \\ast S(t) \n",
    "\\end{equation}\n",
    "\n",
    "#### Exercise\n",
    "Compute the convolution of the source time function 'ricker' with the Green's function of a Vertical displacement due to vertical loads. Plot the resulting displacement."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# call the source time function\n",
    "T = 1/5   # Period\n",
    "src = ricker(dt,T)\n",
    "\n",
    "# Normalize source time function\n",
    "src = src/max(src)\n",
    "\n",
    "# Initialize source time function\n",
    "f = np.zeros(nt)\n",
    "f[0:int(2 * T/dt)] = src\n",
    "\n",
    "#################################################################\n",
    "# COMPUTE THE CONVOLUTION HERE!\n",
    "#################################################################\n",
    "\n",
    "\n",
    "#################################################################\n",
    "# PLOT THE SEISMOGRAMS HERE!\n",
    "#################################################################\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
